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ABSTRACT 
Direct  methods  are  proposed  to  solve  Qv  =  g  where  Q  Is 
a  quasl-trldlagonal  matrix.   Such  matrices  arise  in  ssolvlng 
finite  difference  equations  of  problems  for  pure  elliptic  equations 
and  for  the  symmetric  positive  systems  of  K.  0.  Frledrlchs,   An 
example  is  given  of  the  latter  type  for  solving  a  boundary  problem 
for  the  Trlcomi  equation.   The  methods  employ  a  partitioned 
decomposition  of   Q  into  a  product  of  lower  and  upper  triangular 
matrices  and  a  criterion  is  given  for  the  methods  to  apply. 
Direct  methods  seem  especially  useful  where  iterative  methods  are 
found  not  feasible.   Other  methods  for  special  types  of   Q  are 
proposed  which  Involve  the  inversion  of  relatively  small  matrices. 
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QUASI- TRIDIAGONAL  MATRICES  AND  TYPE- INSENSITIVE 
DIFFERENCE  EQUATIONS 

1.   Introduction. 


In  solving  linear  partial  differential  equations  by  finite 
difference  methods  a  boundary  value  problem  Is  reduced  to  solving 
a  set  of  linear  equations.   In  such  Instances  the  matrix  Involved 
usually  takes  a  special  form  and  consists  mainly  of  zeros.   Many 
of  these  matrices  fall  Into  the  class  to  be  considered  here  which 
may  be  called  quasi- trldlagonal  matrices.   That  is,  we  consider 
partitioned  matrices  of  the  form 


(1.1)  Q  = 


M^  E^  0     .  . 

Dg  M^  E^  0   • 

0   D^  M_^  E,  0 
3   P   3 


0 
0 
0 


^  t^n'  %>  ^n^? 


0   D^  T  M^  ,  E^  , 
q-1   q-1   q-1 


D 


M 


where  the  D  ,  M  ,  E   are  matrices  with  the  same  number  of  rows 
n'   n'   n 

and  E  ,  M  -i  *  D  p  have  the  same  number  of  columns.   We  propose 
to  solve 


(1.2) 


Qv  = 


by  direct  methods. 

Various  discretizations  lead  to  matrices  of  this  type. 
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That  the  usual  finite-difference  approximation  of  certain  boundary 
problems  for  the  Poisson  and  the  bl-harmonic  equation  yield  matrices 
of  the  form  (l.l)  has  been  shown  by  0.  Karlqvist [4 ] ,  A.  F.  Cornock[2], 
L.  H.  Thomas [8]  and  others [6,?].   Whereas  the  usual  method  of 
solution  of  the  resulting  linear  equation  is  by  iteration  these 
authors  propose  direct  methods  for  solving  the  above  problems. 

The  processes  described  here  are,  in  part,  extensions  of  those 
described  by  Karlqvist [4]  and  Cornock[2].   Furthermore  it  is  shown 
that  the  finite  difference  equations  obtained  from  symmetric 
positive  systems,  as  defined  by  K.  0.  Frledrichs[ 5] ,  also  fall 
into  the  class  of  matrices  of  the  form  (l.l).   These  include,  in 
addition  to  pure  elliptic  or  hyperbolic  equations,  a  certain  class 
of  boundary  problems  for  equations  of  mixed  type  such  as  the 
Tricomi  equation.   Where  Iterative  methods  for  such  problems  seem 
to  present  difficulties,  even  when  Q  is  positive  definite  and 
symmetric,  the  direct  methods  for  solving  (1.2)  are  shown  below 
to  be  feasible  for  this  larger  class  of  problems. 

A  criterion  is  given  for  the  process  to  apply  which  is  similar 

to  that  found  for  the  LDU  theorem[3].   It  is  also  shown  that  v/hen 

the   M^   have  the  same  order  p,   and  the   D   are  easily 
n  ^ '  n 

invertible,  then  the  process  may  be  reduced  to  multiplication  of 
matrices  and  the  inversion  of  one  matrix  of  order  p.   This  last 
fact  was  noted  by  Cornock[2]  for  the  Poisson  and  bi-harmonic  case. 
More  generally  if  for   r  =  1,  2,...,  k   the  matrices 


Vl  +  ^^'°"  \'    '^O  =  °'  ^k  =  ^ 
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all  have  the  same  order  p   then  it  is  shown  that  the  process 
may  be  reduced  to  multiplications  and  the  inversion  of  k 
matrices  of  orders  p   ,...,  p,  . 

A  code  to  solve  (1.2)  has  been  written  by  Max  Goldstein 
for  the  l.B.M.-70^  at  N.Y.U.  and  this  code  has  been  successfully 
applied  to  a  symmetric  positive  problem  for  the  Tricomi  equation. 

This  code  has  also  been  applied  to  solving  pure  elliptic 
problems  to  compare  the  direct  method  with  iterative  methods. 
The  direct  method  was  used,  for  instance  to  solve  the  bi-harmonic 
boundary  problem  of  a  simply  supported  rectangular  slab.   A 
comparison  of  running  time  would  indicate  the  direct  method  is 
considerably  faster  than  iterative  methods  for  this  type  of 
problem. 

2.   Direct  methods. 


We  seek  a  reduction  of   Q  to  the  form 

(2.1)  Q  =  LU 

where   L  and   U  are  square  matrices,  partitioned  in  the  same 
manner  as   Q,   of  the  fonn 

(2.2)  L  =  [C^,  I^,  0]J 

(2.3)  U  =  [0,  A^,   E^]^ 


where   I   is  a  unit  matrix  of  the  same  order  as   M,  .   We  also 
n  n 

partition  the  column  vectors  v      and   g   in  the  same  manner. 
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v^ 


S] 

g- 


V  = 


g  = 


Comparing  right  and  left  hand  sides  of  (2.1)  we  set, 
for   i  <  n  <  q, 

A,  =M-,,D   =CA   -,  ,M  =CE  ,    +   k      . 
1  1  '   n    n  n-1  '   n    n  n-1    n 

If  the  A   are  non-singular,  the  A   and  C   may  be 

obtained  recursively  from 


(2.4) 


^n  =  '^n  -  DnA;;-lEr,_i  >    '^i    =  \ 


(2.5) 


n  n  n-1 


C„  =  D„A„",   ,   1  <  n  <  q 


To  solve  for  v   let   Uv  =  y,   then   Ly  =  g;   y   and   v  may  tnen 
be  obtained  recursively  from 


(2.6) 


■^n   ^n 


Vn-1   1  <  n  <  q  , 


where 


Where 


y^^  =  g^  '   a^d 


(2.7)       v„  =  A"^-(y  -E  V  ,,  )   1  <  n  <  q 
^    '       n    n  ^''n  n  n+1 '  2:^        ^  ^ 


V   =  A  y 
q    q  -^q 


We  will  refer  to  this  recursive  method  (2.4)- (2.?)  as  an 
LU- process.   It  should  be  noted  that  this  pi-ocess  I'equii^es  the 
irverslon  of   A  ,  A^ , . . . ,    A   and  their  storage  for  use  on  the 
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"backward  sweep"  (2.7)- 

In  certain  cases  most  of  these  inversions  may  be  avoided. 
For  instance,  if  all  the  matrices   M^,^  are  of  the  same  order  p 
and  the  D   are  easily  invertible,  we  may  premultiply  Q,     by 
the  auasi-dlagonal  matrix 


D, 


-1 


D 


-1 


\ 


/ 


D 


-1 


We  may  thus  assume  that  D^^  =  1,2  <  n  <  q.   If  we  let 

A  =  h""'',  H   then  the  LU-process  becomes 
n    n-j.  n 


(2.8) 


^^n  -  "n-1  ^•''n  "  ^■^n-2  ^n-i   2  <  n  <  q 


where  H^,  =  I  ,  H-j^  =  lA-^    ; 


(2-9)      H^.i(yn-Sn)  =  "  V2yn-1  '   2  <  n  <  o 


where  y^^  =  g^^  and  H^v^  =  \^j,y^    - 

If  we  let  u^  =  ^n-l^n  ^^®"  ^^°^'   (2-9) 


(2.10) 


^n  =  VlSn  -  ^n-l'  ^  <  ^^  ^ 


where  u-,  =  g.  •   To  obtain  the  v   we  go  back  to  the  original 
equation  (1.2)  whence 
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(2-li)     Vl  =  Sn  -  Vn  -  Vn-fl  '   1  <  n  <  q  -  1 

Where   v    =  g   -  m^v   and   v    is  solved  from 
q-i   °q    q  q        q 

(£.18)  W'\- 

We  have,  in  this  case,  only  one  matrix  H   of  order  p 
to  invert. 

This  recursion  (2.8)-(2.12)  will  be  called  an  H-process. 
This  method  may  be  used,  for  instance,  in  solving  problems  for 
the  Poisson  or  bi-harmonic  equations  over  a  rectangle.   In  the 
case  of  Poisson 's  equation  the   q   inversions  of  the  LU-process 
are  reduced  to  solving  only  one  set  of  p  equations  where  p 
is  the  number  of  mesh  points  on  a  horizont-^l  line.   This  fact 
was  noted  for  these  equations  by  Cornock[2]  in  specific  examples. 

For  regions  which  are  made  up  of  rectangles  a  similar  result 
may  be  obtained.   For  instance,  in  solving  Poisson's  equation 
for  an  L-shaped  region  made  up  of  two  rectangles   R, ,  Rp   each 
having  p  ,  p   points  on  a  horizontal  mesh  line,  respectively, 
only  two  matrices  of  order  p,   and   Pp   need  be  inverted. 

To  treat  the  general  case  let  us  assume  that  for 
1  <  Qi  <  qp  <'  •  •<  Qi,  =  Q  >  Qq  =  0   the   M   have  the  same  order 
p   for   qj,_T  <  n  <  q  .   The  matrix  Q  can  then  be  partitioned 
again  by  combining  those   M   having  the  same  order.   That  is,  let 


n'  n'  n  -"l 


where 


-  9  - 


^r-M 


0  ...0  D 


0 


0 


,   1  <  r  <  k  , 


M'   -  =  [D  ,M  ,E  ]  ^1^  ,   0  <  r  <  k 
r+i    ■■  n'  n'  n-'q  +1  '    — 


0 


r+1 


0 

E, 


^^r+1 


0'  •  -c 


0  <  r  <  k  -  1 


The  LU-process  can  be  performed  for  this  new  form  for  Q.. 
Let  g' ,  y '  ,  V,' ,  1  <  r  <  k  be  corresponding  column  vectors  of 
dimension  p   repartitloned  as   Q,   where,  for  instance,  we  denote 


r-h. 


V 


q  +1 
^r 


V 


'r+1 


,  0  <  r  <  k 


Let 


^1 


°1' 


y- 


S2 


C'y'   where 


D- 


P  '  A  ' 
^2^1 


If  we  let 


(2.15) 


y{  =  >'4w| 


w'   may  be  obtained  by  an  H-process  for  (2.1;>).   If  we  denote 
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the   H  matrices  that  enter  here  by   H-,  ,  .  .  .  ,  H    then  only 

H    need  be  Inverted  and  only  w',H   i,  H    and   H 

(or  W   and   H^  H   -,  )   need  be  saved  for  later  use. 
1        q^   q^-1' 


To  obtain  y'   we  must  compute   A' 


not  require  the  Inversion  of 


A' 
^1 


However  this   does 
as   is   seen  from  the  formula 


D'(A|)-^E|    =   D^ 


A 


-1 


A, 


-1 


\*i\l\ 


A 


0" 


-1 


I  Ii 


^i 


I„   is  the  unit  matrix  of  order  p   and   A   ,-,,...,  A 
r  "^r        q  +1      q 


waere 

are  the   A  matrices  that  enter  in  the   LU  decomposition  of 

A'^i  •   Thus   A  •   differs  from   M'   only  in  that   M   ,  -,   is 
r+i  d  d  q-j  +1 

replaced  by 


r+1 


^^+1  -  \+l«q^«q^-lEq^ 


That  is  A^  is  also  of  the  form  (l.l)  and  has  all  of  its  M 
matrices  of  the  same  order  Pp.  We  are  therefore  reduced  to 
the  previous  case.   We  may  proceed  in  this  manner  until  all 
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the  w' ,   1  <  r  <  k  have  been  computed,  where  w' 


k 


V 


k 


To  obtain  the  other  v'   we  note  from  (2.7)  that  for  r  <  k 


V'  =  w 


(A' )'^E'v' , , 
^  r'    r  r+1 


If  we  let 


(2.14) 


A'z' 
r  r 


E'v'  , 
r  r+1 


the   above  -'—process  applies  to  (2.l4)  to  find   z'.   However 
this  H-process  need  not  be  repeated  since 


E'v'-, 
r  r+1 


E  V 
q   q  +1 , 
^r  ^r   i 


That  is,  all  the  vector  components  of  the   "g"  of  this  problem 
except  the  last  are  zero.   It  follows  from  (2=10)  and  (2.12) 
that  to  obtain  z   we  need  only  solve 


W    =  \-l^q/q^+l 


where  H    ,.,..., H   n,H    are  the  H-matrices  corresponding 

Qp.i+i  q^,-!     q. 


to   A'  . 

r 


^r     ^r 
The  other  components   z 


q  1+1' 

^r-l 


,  z    T   are  obtained 


by  multiplications  as  in  (2.11). 

Tiius  we  see  that  for  a  problem  involving  k  sets  of 
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M-matrlces  of  equal  order  only  the  k  matrices  H   ,  H   ,  .  .  .  ,  H 

^1    ^2  ^k 

need  be  inverted.   Thus  In  solving  Polsson's  or  the  bi-harmonic 

equation  over  a  region  made  up  of  k  adjoining  rectangles  the 

inversion  of  at  most  k  relatively  small  matrices  is  required. 

A  similar  statement  applies  in  higher  dimensions. 

We  note  that  in  the  particular  case  where   Q  is  the  discrete 

Laplacian,  where   '^  =  J»  ^   =   1>2,...,  q,  D  ,  E   are  identities 

and   J  =  [1,-4,1]^  ,   the  inverse  of   H   can  be  given 

explicitly [4 ] .   If  we  denote  the  allied  Chebyshev  polynomials  by 


,  /  \    sinh(q+l)x     „    , 

h  (a)  =  ^^ — '—      ,      2   cosh  x  =  a  , 

"^       sinh  X 
then  H  =  h  (j).   The  eigenvalues  of  J  and  H   are  given  by 


\,  =  2  cos|I^-  4,  hq(A^;,  m  =  1,  2,...,  p, 
respectively,  and  the  matrix  of  normalized  eigenvectors  is 

P+^       P+1   k,m=l 
If  we  let 


L  =  [hq{A^),...,  hq(Ap)] 


be  a  diagonal  matrix  then 


H"-"-  =  GL'-'-G. 


For  large  problems  the  H   may  be  ill-conditioned.   For 
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the  above  problem  the  P-condition  P(H  ),   or  ratio  of  largest 
to  smallest  eigenvalue  of   H   is 

TT  sinh  6  (q+l) 
P(Hq)^ , 

(p+l)sinh  6-sinh(^ii  ir) 

so  that  this  may  get  quite  large  for  large  q  =  p=   The 
Inversion  of   H   may  then  present  severe  difficulties.   It  may 
however  be  feasible  to  break  the  problem  into   k  groups  as 
Indicated  above  even  though  all  the  M   have  the  same  order. 

Since  (2.11)  represents  a  marching  process  there  is  also 
the  possibility  of  severe  loss  in  accuracy  in  the  value  of  v-, 
for  large  q.   The  value  of   p   does  not  appear  to  be  an  important 
factor  in  this  loss  for  the  discrete  Laplacian  problem.   In  the 
cases  tried  for  q  =  5,  a  704  code  yielded  results  accurate  to 
at  least  five  significant  digits  for  values  of  p  =  5,  10,  20,  40. 

5.   Criterion  for  decomposition. 

A  sufficient  condition  for  the  validity  of  the  decomposition 

Is  similar  to  that  given  for  the   LDU  theorem[3].   Let 

'ME 
Q   =  M^,  Q   J 

1^2^. 
^k  =  fl^n'^'EjV..,  Qq  =  Q. 

If   Q-,  ,  Op ,  •  .  .  ,  Q   are  non-singular  then  the  decomposition 
(2.1)  -  (2.3)   exists,  and  is  uniquely  given  by  the  recursion 

(2.4),  (2.5).   In  this  case   det  Q  =  TT  det  A   and  if  the   M, 

k=l     ^  ^ 
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all  have  the  same  order  and  D^ , » . » ,  D   are  non-singular  then 
det  Q  =  det  H  , 

The  proof  follows  easily  by  induction.   Prom  the  Schur- 
Probenius  formula 


det  Q-^-i  =  det 


R  M 


q+1 


v-1. 


=  det  Qq  =  det(Mq_^^  -  RQ^  K) 


where  R  =  fO— >0  D  j)  ,   K  = 


0 


0 

L    J 


From  the  Inductive  hypothesis 


Q. 


A' 


-1 


A 


0 


so  that  RQq-'K  =  ^q^iAq'^E^  and  det  Qq_^-j^  =  det  Q^'det  A^^^ , 

This  Implies  that  ^   ,   is  non-singular  and  also  proves  the 
formulas  for  det  Q, 

A  sufficient  condition  for   Q   .  . 
is,  clearly,  that  for  all  u 


Q   to  be  non-singular 


(3^) 


T      T 
yu  Qu  >  u  u 


for  some  non-zero  constant  y. 
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4 .   Applications . 

We  consider,  as  an  application  of  the  above  processes, 
the  problem  of  solving  a  certain  boundary  problem  for  the  Trlcoml 
equation 

(4.1)  Tu  =  yu^^  -  Uyy  -  f (x,y)  =  0 

by  a  difference  approximation  given  by  Frledrichis  [  5]  •   These 
difference  approximations  for  symmetric  positive  systems  have 
been  further  investigated  by  C  K.  ChL.[l]  ^^ 

The  problem  for  (4.1)  is  posed  for  the  parallelogram 

(4.2)  P:  |y-x|  <  t  ,  |xl  <  r  ;  t,  r  >  0 
such  that 

(4.3)  Tu  =  0  ,  (x,y)  e  p 


(4.4)      u^  +  u   =  0  for   |x-y|  =  t  ,  lx|  <  r 


(4.5)      ^v  "  °  ^°^  x  =  -  r  ,  iy-x|  <  t 


No  condition  is  specified  on  x  =  r. 

Since  the  treatment  given  by  Frledrlchs  calls  for  a 
rectangular  region,   P  is  transformed  by   ':=x,  y=y-x 
into  the  rectangle 

(^.6)  Ul  <  r  ,  hi  <  t 

and  the  equation  (4.1)  is  written  as  the  system 
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/y  0^      /y  1  \       'f 
.0  1   J    ^         yi  1    "l      0  / 

Where  y  =  «  +  t],  v  =  (v^  ,v„)  ,  v.,  =  u  ,  v^  =  u,  .    It  is  shown 
in  [i]  and  [5]  that  after  a  premultipiicatlon  of  (4.7)  by  a 
suitable  two  by  two  matrix  of  the  form 


/  p  y 

(4.7)  and  its  boundary  conditions  can  be  brought  into  the 
required  symmetric  positive  form.   That  is  (4.1),  (4.4),  (4.5) 
can  be  written  in  the  form 


(4.8)  |(a^v^  +  rA^  +  (a'v)^  +  (a^v)^)  +  J, 

where  for  a  constant  ^ 

■'  py  ./ 

(4.9)  a^  =         ,  p  =  1  +  e^ 

y  P 

(l^p  ;y  p -y  , 

(4.10)  a^   =  -  ■ 

p  -y   up/ 


V  =  a 


S 


\ 


.pi* 

(4.11)     g  = 


^4.ev  e  \ 
(4.12)     y-  =  -  i(a!  +  a!))  =  ^ 


As  is  shown  in  [l]  and  [5]  the  boundary  conditions  can  be  written 
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In  the  form 


(4.15)  pv  =  ^v 


where 


py  y 

(4.14)     3  =  -  a^,  >a  =  -  a^  =  -(         ^  f or  i  =   i     =  -  r, 

\y    f  -P 


(4.15)  ^  =  a^,  )x  =  aS  =  a^     for   4  =  ^_^  =  r 

(4.16)  P  =  -  c^  ,  ju  =  -  a^  for   11  =  n_  =  -  t  where 

2p^  -  y(p^+l)    (l+p)(p-y) 


-Tl      1 
"   =  P^ 


(l+p)(p-y)       p^  -  2y+l 


(4.17)     3  =  a^  p  =  a^   for  t)  =  n_j_  =  t. 


To  obtain  the  difference  equations  we  divide  the  Intervals 
(-r,r),  (-t,t)  by  the   2p-l,  2q-l   points, 

(4.18)      i^    =    (l-p)A^    1  <  1  <  2p-l,  1  <  P, 


(4.19)      nj  =  (j-q)ATi     1  <  J  <  2q-l,  1  <  q, 


respectively,  where  p  and  q  are  odd  Integers  and 
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A4  =  -— Y  ,    Atj  =  — y  are  the  respective  Interval  lengths. 

According  to  [5]  an  equation  is  written  only  for  the  pq 
odd  numbered  points  {^.,r\    ),    1    =   2m- 1,  J  =  2n-l,  1  <  m  <  p, 
1  <  n  <  q.   If  we  write   v^   =  v(^^,r|  ),  a|  .  =  a^(^^,Ti.)   and 
similarly  for  the  other  variables  then  the  difference  equations 
are  given  for  an  interior  point   1  <  m  <  p,  1  <  n  <  q,   by 


(^■•2°)     k    (4l,J-i.2,J  -  4-l,/i-2,j)  - 


k    (4,j+l^i,j+2  -  «l,j-l^i,J-2)  + 


•^ij^ij    Sij 


where  h  =  2A^,  k  =  2Ari. 

For  a  boundary  point  at  least  one  of  the  subscripts 
1  +  Ij  j  +  1   falls  outside  the  range  prescribed.   In  that 
case  the   a  and   v   of  the  corresponding  term  are  both 
evaluated  at  that  boundary  point  and  then  replaced  according  to 
the  rule  (4.13).   The   2h   and/or  the  2k;   in  the  difference 
involved  are  also  replaced  by   h   and/or  k,   respectively. 

As  an  Illustration  we  consider  the  equation  for  the  point 
i^-\   ''O.Oj   not  a  corner: 


<*•">     H(-l,j",,J-='r/l,) 


■•"  QV  (<^    ^-o-l^l   lO-Q  -  Cl^   .  -,V      „)  +  K.  ,V.    =  g   . 
dK  ±,j  +  l    l,J+2     -i-,J-±    l,J-2'      Ij  l.j    ■=  Ij 


-  19  - 


For  the  corner  pclnt   (^-,  jil-.)   j>  2k   a 


V 


are  replaced  by  1,  k  and  "^^l^ll'   respectively  in  (4.21). 

The  other  types  of  boundary  points  are  treated  In  a  similar  manner. 

Since  each   (m,n)   yields  one  equation  we  obtain  pq   equations 
In  the  pq  unknowns   v.  .,   or  2pq   scalar  equations  for  2pq 
scalar  unknowns.   For  each  of  the  equations  we  note  that   v^ . 
appears  with  at  most  four  of  its  neighbors   v..  ^^  ..  ,  v,  ^  ., 

^^  J-+2  ,  J.     ±-d  ,  ^: 

V.  .   ,-,      V.  .  ^.   If  we  set   v*^  =  v.  .   and  denote  by   v'   the 
i,j-2    i,j+2  m    ij 

vector  with  2p   scalar  components  arising  from  the  points 

(4  ,T]„  _    )      on  a  horizontal  line,  the  difference  equations  take 

the  form 


(^•22)  ^^X"'  -  <  <-.  *^".<*  <  Ci  -  <  C'  -'  ^' 


m 


1  <  m  <  p 


1  <  n  <  q 


This  system  can  be  written  in  the  form  (l.l),  (1-2)  where 


n   ,n   ^n-,p 


^n   =    [^m'^m'V^Ll 


n        "^    '    m'     ■'m=i 


V    = 


E„    = 


n   oiP 


n 


lo.V°lLi 


1~ 


!      1 

1 

1 

1  • 

•      • 

■ 
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and  where  for      l<m<p,    l<n<q 


d^   =  -  i^  a^  ^    .n   _    1^  „^ 


m 


2k   ^l,j-l    '      m   -        2h   ^i-l,j 


%  -  ^ij        '      ^ra  "  2h   °'i+l,J 


^m  -  2k   ^i,j+l      '      ^m        ^IJ 


For   the   boundary  points      n=l,    q,    l<m<p    , 


1    .   o        H^    =   -   1   r^^ 


di    =  0    ,    d^   =  -  T7  a 


m 


m   ~        k    '^i,2q-2 


m       k      1,2  m 


and   for     m=l,    p,      1<"<Q 


n        „  "^   _       i     ^ 

^1    -   "      '      ^p   -   "    h   ^p-2,j 


'1    -  h   °'2j    '      "^p        ^■ 


For     m=l,    p,      l<n<q, 

and  for  corner  points 

m    ^ij    h   ij   k   ij 


J 
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^m  =^1J  -^  h  °^ij   -  k  ""ly      m  =  p,  n  =  1, 


D    =  K  .    .     -     T     O..    .     +    T-     Oi.. 

m    ij   h   ij    k   ij 


m  =  1 ,  n  =  q. 


where  In  all  cases   1  =  2m-l,  J  =  2n-l. 

Thus  this  boundary  problem  for  the  Tricomi  equation  yields 

a  matrix  of  the  form  (l.l).   It  has  been  shown  by  Chu  [l]   that, 

-1        ^ 
if   £,  r,  t   are  properly  chosen,  e.g.,   €  =  /2,  r  =  '^/2 , 

t  =  /5,   this  matrix  satisfies  the  Inequality  (J.l).   The   LU-process 
is  therefore  applicable.   Since  the   M   all  have  the  same  order  we 
may  in  fact  use  the  H-process  providing  the   D   are  non-singular  and 
easily  invertible.   The  D   are  quasi-diagonal,  and  a  simple  compu- 
tation shows  that  for  the  above  choice  of   e,  r,  t,  a^  is  nrn- 
singular . 

The  above  mentioned  code  was  used  to  solve  the  problem  described 
above  for  the  choice  of 

^{^>y)    =   6yx  -  4y^  +  y-l-2x 


and  was  run  for  various  values  of  p  and  q.   The  solution  to  the 
analytic  problem  (^.3)- (^.5)  can  be  given  explicitly  by 


(4.23)     u(x,y)  =  {y-x)Ui  +   x)  -  V25 


£^ 


and  for  (4.8)-(4.17)  it  is  given  by 
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(4.24) 


vi  =  n^  -  n(i  +  20  -  V25 


v„  =  11(1  +  2|) 


For  the  case  of   p  =  13,  q  =  15   the  code  yielded  an 
answer  accurate  generally  to  three  significant  digits.   The  values 
given  at   l(0,0),  II(^,0)   by  tne  code  are  illustrated  by  tables  I, 
II  respectively: 


pxq     J     3x3        5x5        7  X   13        11   X  7      15   X   15 
v^     I     -6.22        -4.07  -4.02        -4.27        -4.04 


I 

-2  ■ 
xlO        " 

^2 

-0.57 

0.08 

0.01 

-0.04 

-0.01 

^1 

-7^13 

-4.96 

-4.34 

-3.99 

-3.99 

II 

_2 
.   xlO 

^2 

-1.01 

0.54 

0.14 

-0.21 

-0.05 

The  exact  solution  given  by  (4.24)  for   t]  =  0   is 


1 


04,  v, 


0, 


Other  problems  for  the  Tricomi  equation  were  run  to  check  the 
influence  of  the  positivity  and  boundary  conditions.   In  one  case 
the  two  by  two  premultipller  matrix,  needed  to  guarantee  "positivity,  " 
was  omitted.   Again  an  approximate  solution  was  obtained  with  a  slight 
loss  in  accuracy.   For  instance  the  value  at   I  given  by  the  code 
was   v^  =  -.0454,  v^  =  -.0019,   for   p  =  11,  q  =  7. 

A  problem  for  the  homogeneous  Tricomi  equation  with  inhomcgeneous 
boundary  conditions   ([j.-p)v  =  f,   for  a  given  f,   also  yielded 
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results  similar  to  those  given  above. 

Symmetric  positive  systems  for  dimensions  higher  than  two 
may  be  treated  in  a  similar  manner. 

As  is  pointed  out  in  [2]  and  [4],  problems  for  the  bi- 
harmonic  equation  can  also  be  put  into  the  form  (l.l).   The 
LU-method  was  carried  out  for  a  simply  supported  rectangular 
plate  and  was  found  to  give  results  accurate  to  about  three 
significant  digits  for  a  30  by  ^0  mesh.   The  running  time  for 
this  problem  was  about  an  hour  on  the  IBM-704. 

The  above  methods  may  also  be  combined  with  iterative  or 
group  relaxation  methods  where  each  individual  group  relaxation 
is  done  by  a  direct  method.   This  has  already  been  proposed  for 
a  multigroup  diffusion  problem  by  Nohel  and  Timlake  [6]. 

It  may  also  be  noted  that,  in  problems  where  higher  order 
difference  schemes  are  available,  the  direct  methods  will  require 
relatively  little  additional  operations  and  thus  one  may  require 
fewer  mesh  points  in  a  given  problem. 
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